On the importance of conserving mass in sea ice models 



m 
o 



Or 

i 

o 

O 

>> 

d' 



> 

(N 
O 

in 
o 

m 



x 

S3 



W. Moon 1 -* and J. S. Wettlaufer 1 ' 2 >t 

1 Yale University, New Haven, CT, 06520, USA 
2 Nordic Institute for Theoretical Physics (NORDITA), 10691 Stockholm, Sweden 

We describe how a long standing approach used in the thermodynamic modeling of sea ice fails to 
conserve mass. The missing mass is traced to a term that is equivalent to neglecting a leading order 
latent heat flux and we demonstrate its influence using energy balance models with a fractional 
ice cover. It is shown that this neglect is particularly acute in a decaying ice cover approaching 
the transitions to seasonal and ice-free conditions. Accordingly, it is suggested that it may be of 
considerable relevance to re-examine the relevant climate model schemes. 



INTRODUCTION 

The original thermodynamic theory coupling sea ice 
and climate dealt with the system as a column of at- 
mosphere, ice and ocean [11]. This approach is the cor- 
nerstone of contemporary theoretical studies [e.g., 5, and 
refs. therein] and it underlies the sea ice components of 
all contemporary climate models [e.g., 3]. We understand 
that the ice cover presents a distribution of ice thicknesses 
g(h) to the atmosphere and ocean that force its growth, 
decay and deformation [20]. Although the treatment of 
this distribution as a continuous diffcrentiablc function 
is based on clear reasoning, its practical implementation 
in either simple or complex models is a major challenge. 

r] developed an implementation scheme for the theory 
of [20] wherein both ice and open water are considered as 
part of a grid cell. In such a scheme, which forms a frame- 
work for the approaches used in a wide range of climate 
models, the areal fraction of ice and water is computed at 
each time step. The scheme emerged at a time when the 
perennial ice state was not questioned. However, it does 
not conserve the total mass in a grid cell and the manner 
in which it violates this conservation law is of particu- 
lar importance as the state of the ice cover changes from 
perennial to seasonal. Here, we demonstrate this in a 
simple model. The question of how, and how rapidly, the 
ice-cover may decay towards the seasonal state, is the 
main implication of the analysis that follows. 



Multiple Sea Ice Cover States 

A main focus of the attempt to discern the origin of the 
decline of the Arctic sea ice cover is the evolution of the 
summer sea ice minimum [e.g., 2, 9] and the associated 
question of whether that minimum shrinks to zero leading 
to a seasonal ice cover in which there is only winter-time 
ice, which melts away each summer. The approaches to 
the problem range from theoretical treatments [1, 5, 6, 
12, 13, 16, 19] and global climate model simulations [e.g., 
8, 21, 22], to syntheses of observations [e.g., 14, 18]. 

The rudiments of the ice-albedo feedback provide the 
framework for examining the nature of transitions from 



the perennial ice state to either a seasonal or ice-free sea 
ice state. In the framework of simplified versions of the 
column model of [11] the ice-albedo feedback treats the 
sea ice albedo as a function of ice thickness h, transition- 
ing continuously from that of sea ice to that of the ocean 
[5, 6, 12]. As the greenhouse gas forcing, modeled as an 
additional surface heat flux AFo, increases these theo- 
retical approaches capture the general conditions under 
which, and the nature of the transitions between, peren- 
nial, seasonal and ice free states. [5] provides a summary 
of the models and methods used to predict four general 
scenarios under which ice retreat may occur as AFq in- 
creases. 



TREATING PARTIAL ICE COVER 
Column Models 

Even in the simplest of models, it is reasonable to at- 
tempt to model partial ice cover, which requires an ocean 
mixed layer that is always in communication with the at- 
mosphere unless the ocean is completely ice-covered. The 
most common methodology to evolve ice area A appears 
to have originated from that of [7] . This approach is the 
same regardless of the degree to which a model handles 
g(h) and the mechanical and thermodynamical contri- 
butions to it. Thus, it is prudent to take the minimalist 
approach to demonstrate the matters at hand. A detailed 
derivation of a simple ice column model from that of [11] 
is given by [6] and we summarize the relevant aspects 
here. 

When the temperature of the ice Tj < °C it evolves 
along with the ice thickness according to 
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where L is the latent heat of fusion, c is the effective 
heat capacity, and k is the thermal conductivity. The 
sum of sensible, latent, downward and upward longwave, 
and shortwave heat fluxes at the surface is given as F top , 
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and all but the upward longwave flux are specified from 
observed radiation climatology. The flux from the ocean 
mixed layer into the base of the ice is Fb- 

When Ti = °C the temperature evolves along with 
the ice thickness according to 



L 



dt 
dh 



0. 
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top 



(3) 
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The ocean mixed layer is treated as a thermodynamic 
reservoir with an observationally based characteristic 
depth of H m i = 50m. The flux of heat entrained through 
the bottom of the mixed layer is specified to be F ent = 0.5 
W m~ 2 based on observations. The turbulent heat 
flux between the ocean and the underside of the ice is 
given by Fb = pCp W ChU+ AT, in which p and c pw are 
the density and specific heat of seawater, Ch = 0.06 is 
the heat transfer coefficient, AT = T m i, and u* is the 
square root of kinematic stress at ice-ocean interface, also 
known as the friction velocity. Using a typical observa- 
tional value of u+ = 0.5 cm s" 1 leads to Fb = 7T m j with 
7 = pc pw ChU*o = 120Wm~ 2 . 

Therefore, according to whether Tj < 0°C, a contin- 
uously evolving seasonal cycle is captured, using (l)-(2) 
or (3)-(4). 



Modeling Partial Sea Ice Cover 

Now we proceed to the matter of modeling partial ice 
cover through an extension of column models. First, we 
summarize the approach to this problem developed by [7] . 
Second, although such an approach can be rationalized 
for the perennial ice cover for which it was developed, 
we discuss its quantitative acuity when the ice fraction 
decreases, such as when dealing with climate transitions 
in the overall state of the ice cover. 

7] introduced a methodology to evolve ice concentra- 
tion A, or the fraction of a grid box covered by ice. This 
requires a form of homogenization over the subgrid-scale 
to account for open water. As noted in the introduc- 
tion, this method provides the framework for parameter- 
izations of subgrid-scale open water and thickness dis- 
tributions in many contemporary GCMs with the most 
sophisticated sea ice representations. For clarity we dis- 
cuss the approach in terms of a model that includes a 
single grid box. 

The ice concentration increases when T m \ reaches zero 
and continues to cool so that the mixed layer flux imbal- 



dA _ F ni 
dt Lh 



(5) 



ascribing volume to it. Thus, area increases only when 
the mixed layer freezes, but once it does so that new vol- 
ume of ice only increases by increasing the ice thickness 
at fixed area. Because, with the framework of column 
models, sea ice growth rate is calculated (or specified) 
as a function of ice thickness and season, the value of 
ho controls the rate at which ice cover grows. Although 
the growth rate in winter decreases by a factor of four 
as open water solidifies to a thickness of 50 cm, the ice 
concentration in Hibler's approach increases based on the 
growth rate for open water ~ 12 cm day" 1 (sec figure 3 
of [7]). The value of ho used in [7] is 50 cm. Importantly, 
in this and similar models, the open water fraction is not 
meant to represent an entirely ice-free region. Rather, 
the model domain is split into a fraction containing thick 
ice, with the rest covered by a mixture of open water 
and thin ice, such as in leads. The volume of this thin 
ice is assumed to be negligible compared to the thick ice 
volume, which as we shall see is one of the problems in 
dealing quantitatively with processes such as ice-albedo 
feedback. 

According to energy balance, area decays in this model 
when volume ablates. Hence, when tt- < area de- 
creases as 



dA _ A_dV 
~dt ~ 2V~dt' 



(6) 



The proportionality between volume and area rates of 
change is based on an argument about the ice thickness 
distribution in the model domain as follows [7]. Assume 
that (a) the ice is linearly distributed in thickness be- 
tween and 2V/A, thereby giving a mean thickness of 
V/A, and (b) all of this ice melts at the same rate. As il- 
lustrated schematically in Fig. 1, this gives a rate of area 
decay as the rate of thickness decay times the inverse 
slope of the thickness distribution; 



a , dA 
AA = Ah — 
dh 
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(7) 



An "equivalent thickness" ho is assigned to the new area 



We note here that ice growth is a nonlinear function of 
thickness and here it is computed under the assumption 
that all ice within A is of the mean thickness V/A as 
opposed to the linear distribution between and 2V/A 
used in ablation. 

Finally, due to persistent convergence and divergence 
of the wind field and the net average observed annual 
export of 10% of the ice area, the ice dynamics are rep- 
resented in such a model by requiring that A < 0.95 
and a term —voA to the area evolution equation (which 
accounts for volume export naturally in when mass con- 
servation is obeyed). Export is included in the results 
shown figure 2, but to avoid the clutter in the theoreti- 
cal development we omit the term in the equations that 
follow. This has no effect on the main points. 

Using such a scheme one can derive a partial ice cover 
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model from the column treatment of § as follows. We 
evolve ice volume rather than ice thickness. In the ice- 
covered fraction of the model domain A, this vertical 
thermodynamic growth of the ice is represented by re- 
writing (2) and (4) as 



dt V h 



(8) 



and 



L^=A(-F top -F B ). (9) 

The total heat flux into the mixed layer is written as 

F ml = - AjTmi + F ent (10) 
+ (1 - A) \-F lw (T ml ) + (1 - a m i)F sw ] , (11) 

where F[ w (T m i) is the net surface longwave radiation bal- 
ance which depends on the mixed layer temperature T m i , 
the shortwave radiative balance is F sw , and the albedo of 
the mixed layer is a m i [4, 6]. Therefore, if T„a > 0, this 
leads to heating or cooling according to 



dT m i 
at 



(12) 



and no new ice area is formed, F n i = 0. However, when 
the mixed layer reaches the freezing temperature (T m ; = 
0), supercooling is prohibited such that dT m i/dt = 0, 
and any additional heat loss is available to form new ice; 

{.Fni — Fml*). 



ice cover when the mass of thin ice is small relative to the 
thick ice, it becomes questionable as the area of thin ice 
increases in magnitude and the rate of change of its area 
becomes significant. Indeed, as the ice cover transitions 
to a seasonal state undergoing rapid and large changes in 
areal coverage, and during a part of such a seasonal cy- 
cle when A becomes vanishingly small, the concept of an 
average thickness defined as h = h/A, exhibits obvious 
pathologies. 

Requiring mass conservation modifies the left hand 
side of Eqs. 8 and 9 through the addition of a term 



VdA 
'A dt 



Now, to assess the importance of such a term, 
it is prudent to render these equations dimensionlcss 
through the introduction of the following scalings; 



V = Wo, A = AA , 
L 



t = Tr, 



S 



CpAT 1 



Ti = %AT, 
AT 

F c = k—, (13) 

"0 



is the threshold thickness mentioned 

Ao 



where ho 

above, c p is the specific heat capacity of the ice, AT 
is the temperature difference over ice of thickness ho and 
hence F c is the associated conductive heat flux. The di- 
mensionless ratio S is a Stefan number, which represents 
the relative importance of latent heat to the specific heat 
in the ice, and is large (> 10) here. These scalings lead 
to dimensionless versions of Eqs. 8 and 9 modified for 
proper mass conservation as 
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dV 
df 



V_dA 
Adf 



-%^-F B Y (14) 



Conservation and Nonconservation of Mass 

The principal issue here is that Eqs. 8 and 9 do not 
conserve total mass and hence the evolution of the areal 
fraction of ice is incorrect. The nub of the matter is the 
appropriate grid homogenization of the mixture theory. 
We write the volume of ice in a grid cell of area Aq as 
hAc and use the logic of [7], that there are two ice cat- 
egories, thick (hi) and thin (ho), the latter intended to 
represent both open water and ice up to a value of ho- 
Then, one would naturally write the average ice thickness 
as h = hiA+ho(l — A) and hence the growth rate in a grid 
cell would be h = hiA+ho(l — A) + (hi — ho) A, where the 
dot's denote differentiation with respect to time. How- 
ever, the scheme described in § implements the assump- 
tion that the mass of thin ice is small relative to that of 
thick ice in the following manner. On the one hand, this 
assumption imposes ho — in order to compute the aver- 
age "thick ice" thickness as h = h/A. On the other hand, 
the reality that thin ice grows rapidly is implemented by 
computing the growth rate as ft. = hA+ho(l — A) with the 
term (hi — ho) A set to zero (see Eq. 15 of [7]). Although 
such an approach has merits in the central perennial pack 



and 



dV 
df 



V_dA 
AdT 



A(-Ftop-FB), (15) 



where the fluxes J- are just the dimensional fluxes scaled 
by F c . The importance of developing a non-dimensional 
set of equations for the freezing and melting process is 
to facilitate simple analysis and interpretation of general 
features. Because cS > 1 then the balance of both Eqs. 
14 and 15 require that 



dV 
df 



«1, 



V dA~ 
Adf_ 

dV dA 
* V ~ ~A ' 



(16) 



showing that the neglected term is of the same order 
of magnitude of that kept in the scheme. Importantly, 
in transitioning to a seasonal ice state driven by the ice- 
albedo feedback, correctly capturing the rapidly changing 
areal evolution is crucial. The unphysical nature of 
the methodology which leads to Eqs. 8 and 9 manifests 
itself in untoward parameter regimes in which ice loss 
is expected. For example, when [4] neglected this term, 
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he found that in order to obtain multiple sea ice states 
under greenhouse gas forcing of 1.5 times the present 
value he had to artificially decrease the latent heat of 
fusion of ice by factors ranging from four to ten. Thus, 
for a given radiation balance, this reflects the artificially 
large latent heat flux associated with the neglect of the 
extra term in Eqs. 14 and 15. Moreover, when using the 
correct thermophysical constants, he had to use 5 times 
(16 times) the present greenhouse gas concentration to 
transition from perennial to seasonal ice (multiple sea 
ice states). 



radiative fluxes as 

AF t = [ A [l - a(h)]F s (t)dA' - A[l - a(h)}F s (t) 
Jo 



-F s (t)x 



!-^ln 
2h 



cosh 



tanh 



h a 



(18) 
(19) 



Therefore, the lateral oceanic heat flux ablating the ice 
is h^rr — AFi. Finally, in order to conserve energy flux 
balance at each time step during the melt season, we 
subtract the lateral oceanic heat flux from the evolution 
of ocean sensible heat, thereby avoiding an anomalous 
increase in ocean heat content. 



Energy flux conservation 



During the melt season the contribution to the volume 
evolution of the lateral melting term is insufficient 
to conserve energy flux. In Hibler's scheme, this lateral 
melting is calculated indirectly in order to maintain the 
functional form of the model sea ice thickness distribu- 
tion, save for the constant difference owing to the change 
of the mean. Because greater than 90% of the incident 
sunlight is absorbed by open water, we understand that 
the partitioning of the ablation of the ice cover between 
top, bottom and lateral boundaries is, among other fac- 
tors, a complex function of the open water fraction and 
ice floe perimeter [e.g., 10, 15]. The impasse faced by 
the Hibler [7] model is discussed by Steele [15] in terms 
of the lack of an explicit equation for ice floe perime- 
ter. However, it is still natural to ask for the origin of 
the heat source for lateral melting within this commonly 
used model framework. 

We seek the partitioning of vertical and lateral oceanic 
heat fluxes to account for the contribution h^rr. When 
the average thickness h is used to determine the heat flux 
required to balance the volume change originating in ver- 
tical ablation A^r, the result differs from the analogous 
procedure in which the ice thickness is distributed evenly 
from to 2h. Thus, we conclude that part of this heat 
flux difference AFi over an ice area A is used for lateral 
melting, and is written as 

AFi = / Fi(h)dA' - AFi(h), (17) 
Jo 

where Fi(h) is the net heat flux over sea ice of thickness 
h and A' = — ^h + 2h derives from the sea ice thickness 
distribution. For simplicity of exposition, we view the 
principal contribution to AFi as arising from shortwave 



DISCUSSION 



Now that the essential point has been made using this 
simple analysis, in Figure 2 we show the dramatic effect 
of employing a scheme that does not obey mass conser- 
vation. The missing-mass-term, as described in the ar- 
gument leading to equation 16, has the basic effect of 
neglecting a leading order latent heat flux. Thus, under 
the same radiative forcing, the consequences of missing 
latent heat flux are clearly laid bare. First, because the 
effective latent heat flux is smaller, the volume of ice in 
steady state with the same radiative balance can be larger 
when mass is not conserved as seen in clearly in the top 
panel. The most distinct case is for AFq = 0, where the 
maximum sea ice thickness is 3.5m (2 m) when mass is 
not conserved (mass is conserved). Second, when mass 
is not conserved a larger value of AFq (~5W m~ 2 ) is 
required before both the perennial and seasonal ice states 
vanish. Third, the range of AFo over which seasonal sea 
ice exists in a stable state is practically non-existent (~ 5 
W m~ 2 ) when mass is not conserved (mass is conserved). 
Finally, since the ice cover vanishes at smaller values of 
greenhouse forcing when mass is conserved the heat con- 
tent of the exposed mixed layer is accordingly larger. We 
note that in the original treatment of Hibler [7] there was 
no mixed layer; the ocean temperature was constrained 
to lie on the freezing point and any excess heat absorbed 
was immediately added to a basal heat flux and applied 
to the underside of the thick ice. This treatment is clearly 
unrealistic in the limit of a vanishing ice cover, but the 
missing term in such a scenario becomes all the more im- 
portant, and a direct application of this scheme simply 
amplifies the differences shown in figure 2 so we do not 
include these figures here. 
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FIG. 1. Schematic illustrating the proportionality between 
the rate of change of ice area A and the thermodynamic de- 
crease of volume following [7] . 



CONCLUSION 

We describe how a long standing approach used in 
the thermodynamic modeling of sea ice fails to conserve 
mass. By deriving energy balance models with and with- 
out mass conservation we demonstrate the sensitivity of 
the results to the missing mass and find that in a decay- 
ing ice cover approaching seasonally ice free conditions 
is particularly poorly predicted. Although we have not 
independently analyzed how this erroneous treatment of 
mass conservation has propagated through the range of 
GCMs used, our analysis indicates the possibility that it 
could in fact be one of the underlying features responsible 
for the observed recent Arctic sea ice decline being more 
rapid than is forecasted [e.g., 17]. Thus, it is suggested 
that it may be of considerable relevance to re-examine 
climate model schemes to revise them accordingly. 
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